Phase transitions in self-gravitating systems. 
Self-gravitating fermions and hard spheres models 



Pierre-Henri Chavanis 

Laboratoire de Physique Quantique, Universite Paul Sabatier, 
118 route de Narbonne 31062 Toulouse, France. 

Abstract 

We discuss the nature of phase transitions in self-gravitating systems both in the mi- 
crocanonical and in the canonical ensemble. We avoid the divergence of the gravitational 
potential at short distances by considering the case of self-gravitating fermions and hard 
spheres models. Depending on the values of the parameters, three kinds of phase transi- 
tions (of zeroth, first and second order) are evidenced. They separate a "gaseous" phase 
with a smoothly varying distribution of matter from a "condensed" phase with a core-halo 
structure. We propose a simple analytical model to describe these phase transitions. We 
determine the value of energy (in the microcanonical ensemble) and temperature (in the 
canonical ensemble) at the transition point and we study their dependance with the de- 
generacy parameter (for fermions) or with the size of the particles (for a hard spheres gas). 
Scaling laws are obtained analytically in the asymptotic limit of a small short distance 
cut-off. Our analytical model captures the essential physics of the problem and compares 
remarkably well with the full numerical solutions. We also stress some analogies with 
the liquid/gas transition and with the Blume-Emery-Griffiths (BEG) model with infinite 
range interactions. In particular, our system presents two tricritical points at which the 
transition passes from first to second order. 

1 Introduction 

The statistical mechanics of self-gravitating systems turns out to be very different from that of 
other, more familiar, many-body systems like neutral gases and plasmas due to the unshielded, 
long-range nature of the gravitational force [JJ. Because of this fundamental difference, the 
notion of equilibrium is not always well-defined and these systems exhibit a non trivial behavior 
with the occurence of phase transitions associated with gravitational collapse. If the particles 
are treated as classical point-masses, it can be shown that no global entropy maximum exists, 
even if the system is restricted within a box so as to prevent evaporation 0, [J. A self-gravitating 
system can increase entropy without bound by developing a dense and hot "core" surrounded 
by a dilute "halo" . There exists, however, local entropy maxima (metastable equilibrium states) 
if the condition A = —-§^z < 0.335 is satisfied, i.e. if the energy E is sufficiently large (for 
a given box radius R) or if the radius is sufficiently small (for a given energy E). Since these 
equilibrium states are only local entropy maxima, the question naturally emerges whether they 
are long-lived or if they will collapse to a configuration with higher entropy. In any case, a 
phase transition must occur for A > A c = 0.335 since the entropy has no extremum at all above 
this threshold ||. In that case, the system is expected to collapse indefinitly towards a state of 



higher and higher central concentration and temperature. This is the celebrated "gravothermal 
catastrophe" ||. 

However, if we introduce a repulsive potential at short distances, complete core collapse is 
prevented and it can be proved that a global entropy maximum now exists for all accessible 
values of energy. This effective repulsion can be introduced in many different ways but the 
physical results are rather insensitive to the precise form of the regularization. For example, 
we can study the case of self-gravitating fermions for which an exclusion principle imposes an 
upper bound on the distribution function ||, |f| || 0|. Alternatively, we can consider a classical 
hard spheres gas by introducing an "excluded volume" around each particle || [| . Other forms 
of regularization are possible [It], |Tl|] . When such regularizations are introduced, it is possible 
to evidence properly three kinds of phase transitions of zeroth, first and second order. They 
separate a "gaseous" phase, which is independant on the small scale cut-off, from a "condensed" 
phase in which the particles are tightly bound. This is similar to the liquid/gas transition in an 
ordinary fluid. However, for long-range systems such as self-gravitating systems, the statistical 
ensembles are not interchangeable and phase transitions can occur both in the canonical and 
in the microcanonical ensemble. This results in the existence of two tricritical points, one in 
each ensemble. In that respect the self-gravitating Fermi gas shares some analogies with the 
Blume-Emery-Griffiths (BEG) model with infinite range interactions |T2[ . 

The object of this paper is to provide a detailed description of phase transitions in self- 
gravitating systems. In the first part of the paper, we consider the case of self-gravitating 
fermions. The equilibrium phase diagram was calculated in an earlier paper || and we complete 
this study by determining explicitly the values of energy (in the microcanonical ensemble) and 
temperature (in the canonical ensemble) at which the phase transitions occur. We also propose a 
simple analytical model to describe these phase transitions. The "gaseous" phase is modeled by 
a classical homogeneous sphere while the "condensed" phase is made of a completely degenerate 
nucleus surrounded by a hot atmosphere with uniform density (restrained by the box). The 
mass M* of the nucleus is determined by maximizing the entropy (resp. free energy) versus M* 
for a given total mass and energy (resp. temperature) of the configuration. Quite remarkably, 
this simple model can reproduce the main features of the numerical study. It also allows us to 
determine analytically how the energy or the temperature at the transition points depend on the 
degeneracy parameter. In the second part of the paper, we extend our analytical model to the 
case of a classical gas with a short distance cut-off. This model has been studied numerically by 
Aronson & Hansen || and Stahl et al. || , and our analytical model gives a good agreement with 
their numerical results. It is also consistent with the toy models of Lynden-Bell & Lynden-Bell 



13| and Padmanabhan Jy. The possible astrophysical applications of our study are discussed 
in Sec. IETd*. 



2 Statistical mechanics of self-gravitating fermions 
2.1 The Fermi-Dirac distribution 

We consider a system of N fermions interacting via Newtonian gravity. These particles can 



be electrons in white dwarf stars [H], neutrons in neutron stars |T^, |j, massive neutrinos in 
dark matter models [16, P etc... We assume that the mass of the configuration is sufficiently 
small so as to ignore relativistic effects. Let f(r, v,i) denote the distribution function of the 
system, i.e. /(r, v, t)d 3 rd 3 v gives the mass of particles whose position and velocity are in the 
cell (r, v; r + <i 3 r, v + <i 3 v) at time t. The integral of / over the velocity determines the spatial 
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density 



(1) P = J fd 3 v, 
and the total mass of the configuration is given by 

(2) M = I pd 3 r, 



where the integral extends over the entire domain. On the other hand, in the meanfield ap- 
proximation, the total energy of the system can be expressed as 

(3) E = I f fv 2 dhd\ + lf p$d 3 r = K + W, 



2 J 2 . 

where K is the kinetic energy and W the potential energy. The gravitational potential $ is 
related to the star density by the Newton-Poisson equation 

(4) A$ = AnGp. 
Finally, the Fermi-Dirac entropy is given by the formula 

(5) S = - fit ln^ + ( 1 - iJ i _ L) \dhd 3 v, 

J {Vo Vo V VoJ V VoJ J 

which can be obtained by a standard combinatorial analysis. In this expression, rj Q is the 
maximum value accessible to the distribution function. If g = 2s+l denotes the spin multiplicity 
of the quantum states, m the mass of the particles and h the Planck constant, one has by virtue 
of the Pauli exclusion principle i] = gm 4 /h 3 . An entropy of the form (J|) was also introduced by 
Lynden-Bell [|T7| in the context of collisionless self-gravitating systems (e.g., elliptical galaxies, 



dark matter) undergoing a "violent relaxation" by phase mixing |1| [19], |20[. In that context, 
rjo represents the maximum value of the initial distribution function and the actual distribution 
function (coarse-grained) must always satisfy / < tjq by virtue of the Liouville theorem. This 
is the origin of the "effective" exclusion principle in Lynden-Bell's theory which has obviously 
nothing to do with quantum mechanics. In reality, the mixing entropy introduced by Lynden- 
Bell is a complicated sum of Fermi-Dirac entropies for each phase level constituting the initial 
condition. For simplicity, we shall restrict ourselves to the single level approximation for which 
the mixing entropy coincides with expression (jj|). 

At statistical equilibrium, the system is expected to maximize the Fermi-Dirac entropy at 
fixed mass and energy. Introducing Lagrange multipliers to satisfy these constraints, we find 
that the critical points of entropy correspond to the Fermi-Dirac distribution 

Vo 



(6) / 



l + Ae ^+*) : 



where A is a strictly positive constant insuring that f < Vo an d (3 is the inverse temperature. 
In the fully degenerate limit / ~ 770, this distribution function has been extensively studied in 
the context of white dwarf stars in which gravity is balanced by the pressure of a degenerate 
electron gas [ I4"| . In the non degenerate limit / <C 770, the Fermi-Dirac distribution reduces to 



the Maxwell-Boltzmann distribution 

(7) / ~ ^e-*4+*), 

A 



3 



so that we expect to recover the properties of classical isothermal spheres at low densities H . In 
particular, the Fermi-Dirac spheres, like the isothermal spheres, have an infinite mass and one 
is forced to confine the system within a box of radius R. Physically, this confinement is justified 
by the realization that the relaxation is incomplete so that the conditions of applicability of 
statistical mechanics, regarding for example the ergodic hypothesis, can be fulfilled only in a 
limited region of space. In addition, an astrophysical system is never completely isolated and 
R could represent the typical radius at which the system interacts with its neighbors. 

Kinetic equations have been proposed to describe the relaxation of self-gravitating systems 
towards the Fermi-Dirac distribution ETJ]. If we assume that the system is subject to tidal 
forces and if we allow high energy particles to escape the system when they reach an energy 
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e = \ + $ > e m , an extension of the Michie-King model taking into account the degeneracy 



can be deduced from these equations [21]. For e < e m , one has 



(8) f = Vo 



A + e-P £ 



while / = for e > e m since the stars have been removed by the tidal field. When A — > +oo, 
we recover the Michie-King model E| and when e m — ► +oo, we recover the Fermi-Dirac 



distribution (g). The density associated with this distribution function goes to zero at a finite 
radius, identified as the tidal radius. Therefore, the configuration has a finite mass. This 
distribution function could describe elliptical galaxies and galactic halos limited in extension 
as a consequence of tidal interactions with other systems (I7|, [RJ . This model is of course more 
realistic than the box model. However, in order to exhibit phase transitions in self-gravitating 
systems, the box model provides a more convenient theoretical framework and we shall use it 
in the sequel. 



2.2 Thermodynamical parameters 

The thermodynamical parameters for Fermi-Dirac spheres in the meanfield approximation 
have been calculated by Chavanis Sz Sommeria |J and we shall directly use their results. 
The equation determining the gravitational potential at equilibrium is obtained by substi- 
tuting the Fermi-Dirac distribution @ in the Poisson equation (f|), using Eq. ([]]). Intro- 
ducing the variables ip = /3($ — $q), where $o is the central potential, k = Xe^° and 
£ = (167r 2 V2Gr] /f3 1/2 ) 1/2 r, it can be written 

where I\/2 denotes the Fermi integral 

r+oo n 

(10) /„(*) = / ^~—dx, 



1 +te x 

of order n = 1/2. The boundary conditions at the origin are 

(11) ^(0) = ^(0) = 0. 

In the case of bounded spheres, one must stop the integration of Eq. (Q) at £ = a with 

(12) a --(^^)'\. 
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The parameter a is related to the temperature and to the energy by the equations 
(13) r] = — — = aip k (a), 



Moreover, a and are related to each other by the relation 

(15) aV fc (o) = /i 2 , 
where /i is the "degeneracy parameter" 

(16) ll = ri V5l2ir 4 G 3 MR 3 . 

For a given value of \x and fc, we can solve the ordinary differential equation (P) until the value 



^ = a at which the condition (T5|) is satisfied. Then, Eqs. (|13[)(|14l) determine the temperature 
and the energy of the configuration. By varying the parameter k (for a fixed value of the 
degeneracy parameter lm), we can cover the whole bifurcation diagram in parameter space ||. 
The entropy of each configuration is given by (see Appendix A) 



(17) ^ = - ? -Ar] + ifc(a) + v + In k - ^I s/a (*e*W; 



and the free energy by 

(18) J = S-/3E. 



Note that Eq. ( jlq) is the free energy F = E — TS up to a negative proportionality factor. In 
the microcanonical ensemble, a solution is stable if it corresponds to a maximum of entropy S 
at fixed mass and energy. In the canonical ensemble, the condition of stability requires that 
the solution be a maximum of free energy J at fixed mass and temperature. It can be shown 
0] that this meanfield approach is exact in a thermodynamical limit such that — > +oo with 
/i, 77, A fixed. This implies in particular that N 1 ^ 3 R, TN~ 4 ^ 3 , EN~ 7 ^ 3 and SN^ 1 approach 
a constant value for A^ — > +00. The usual thermodynamical limit N,R —> +00 with N/R 3 
constant is of course not relevant for non extensive systems. 

3 Phase transitions in a self-gravitating Fermi gas 

3.1 The non degenerate limit (// = 00) 

Before considering the case of an arbitrary degree of degeneracy, it may be useful to discuss 
first the non degenerate limit corresponding to a classical isothermal gas (h — > 0). In that case, 
the thermodynamical parameters are given by 

(19) 7] = onj)'{at), 

3 1 e -^(") 

(20) A= 2^(^)~^)2' 



5 



S 1 

(21) — = --lnr/ - 2 In a + ip(a) + rj - 2Ar), 

where a = (AnGfipoY^R * s the normalized box radius and ip is the normalized gravitational 



potential solution of the Emden equation [14 



with boundary conditions 

(23) ^(0) = f(0) = 0. 

The equilibrium phase diagram (E, T) is represented in Fig. [L| The curve is parametrized by 
a which can be considered as a measure of the central concentration. An equivalent parametriza- 
tion is provided by the density contrast 1Z = p(0)/p(R) (see Fig. Q). m * ne microcanonical 
ensemble, the solutions on the upper branch of Fig. |l] (until point MCE) are stable and corre- 
spond to local entropy maxima. The rest of the spiral corresponds to unstable saddle points. 
For A > A c = 0.335, there are no critical points of entropy and the system collapses indefinitely. 
This is the so-called "gravothermal catastrophe" . In the canonical ensemble, only the solutions 
prior to the point CE are stable. They correspond to local maxima of free energy. The rest 
of the spiral corresponds to unstable saddle points. For rj > rj c = 2.52, there is no hydrostatic 
equilibrium and the system undergoes an "isothermal collapse" . These stability results can be 
deduced from the turning point analysis of Katz [23| or by solving the eigenvalue equations 



associated with the second order variations of entropy or free energy p^j . For A < A c and 
i] < rj c , the stable solutions are only metastable. There are no global maxima of entropy or free 
energy for classical point masses in gravitational interaction M. We also note that the region of 
negative specific heats (between points CE and MCE) is stable in the microcanonical ensemble 
and unstable in the canonical ensemble, where it is replaced by a phase transition (an "isother- 
mal collapse" ) . This is expected on physical grounds since it can be shown quite generally that 
the specific heat must be positive in the canonical ensemble |TJ. These results clearly indicate 
that the statistical ensembles are not interchangeable in the case of self-gravitating systems 
(see, e.g., Ref. |26[), contrary to normal matter in which the energy is an extensive parameter. 



In Fig. |3|, we have ploted the energy and the entropy as a function of the central con- 
centration a. We observe that the peaks occur for the same values of a. This is also clear 
from Fig. |] where we have represented the entropy as a function of the energy. When several 
critical points of entropy exist for the same energy, only the one with the largest entropy is 
an entropy maximum. The other critical points are unstable saddle points. Therefore, if the 
system is initially prepared on a saddle point, we expect a transition to occur from a state 
of low entropy to a state of higher entropy. This is not really a phase transition but just an 
instability. A similar diagram has been found for isothermal spheres described in the context 



of general relativity ]2"7 |. In this analogy, the mass-energy M plays the role of the classical 
energy E and the binding energy Eu n d. = (M — Mo)c 2 , where Mq = Nm is the rest mass and 
N the baryon number, plays the role of the classical entropy 5*. In Figs. [|-|6|, we have ploted 
the corresponding diagrams in the canonical ensemble. Here again, the peaks of temperature 
correspond to the peaks of free energy. The free energy has an additional peak which has no 
counterpart in the temperature diagram. However, this peak is not associated to an instability 
and the interpretation of the curves is the same as in the microcanonical ensemble. 
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Figure 1: Equilibrium phase diagram for classical isothermal spheres. For A > A c or rj > r) c , 
there is no hydrostatic equilibrium and the system undergoes a gravitational collapse. 




Figure 2: Density contrast of classical isothermal spheres as a function of energy. The series 
of equilibrium becomes unstable in the canonical ensemble for 1Z > 32.1 and in the canonical 
ensemble for 1Z > 709. The value of energy at which the density contrast tends to +oo is 
A = 1/4 (singular sphere). 
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Figure 3: Entropy and energy as a function of the central concentration a. The peaks of energy 
and entropy occur for the same values of a. 
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Figure 4: Entropy versus energy for classical isothermal spheres. When several solutions exist 
for the same energy, the states with low entropy are unstable saddle points. They can either 
evolve towards the metastable state with highest entropy (see arrow) or collapse to a state of 
ever increasing entropy, as suggested in Ref . [FZUJ . 
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Figure 5: Free energy and inverse temperature as a function of the central concentration a. 
The peaks occur for the same values of a. The free energy presents an additional extremum 
but it is not associated to an instability. 
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Figure 6: Free energy versus inverse temperature for classical isothermal spheres. The inter- 
pretation is the same as in Fig. [|. 
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3.2 Large values of the degeneracy parameter (fi = 10 5 ) 

We now consider the case of self-gravitating fermions characterized by a degeneracy parameter 
fi. We first discuss the case of large values of the degeneracy parameter. The extension of the 
classical diagram of Fig. |l] is reported in Fig. [7]. We see that the inclusion of degeneracy has 
the effect of unwinding the spiral. The evolution of the density contrast along the series of 
equilibrium is depicted in Fig. || In the range A*(/i) < A < A c , there exists several critical 
points of entropy for each single value of energy. The solutions on the upper branch of Fig. 
[7| (points A) are non degenerate and have a smooth density profile; they form the "gaseous" 
phase. The solutions on the lower branch (points C) have a "core-halo" structure with a massive 
degenerate nucleus and a dilute atmosphere; they form the "condensed" phase. According 
to the criterion of Katz [p3fl , specifically discussed in Ref. @ for self-gravitating fermions, 
these solutions are both entropy maxima (EM) while the intermediate solutions (points B) are 
unstable saddle points (SP). These points are similar to points A, except that they contain a 
small embryonic nucleus (with small mass and energy) which plays the role of a "germ" in the 
langage of phase transition. The density profiles of these solutions are given in Ref. ||. 

To be more precise, we have ploted the entropy of these solutions as a function of energy 
in Fig. |9|. The entropy of the unstable phase (points B) is always smaller than the entropy of 
the stable phases, as it should. There is now a crossing point in the diagram, at A = A t (/J,), 
which marks the onset of a phase transition. At that point, the "gaseous" phase and the 
"condensed" phase have the same entropy. As A is increased across the transition point, the 
non degenerate solutions (points A) pass from global to local entropy maxima. Inversely, 
the degenerate solutions (points C) pass from local to global entropy maxima. We expect 
therefore a phase transition to occur from the "gaseous" phase to the "condensed" phase when 
A = Af(fi). The "kink" in the curve S(E) at the transition point where the two branches 
intersect corresponds to a discontinuity of temperature in the equilibrium phase diagram (see 
the vertical plateau in Fig. ^). The specific heat is also discontinuous at that point and turns 
from positive to negative. According to Stahl et al. |J, this phase transition could be called a 
"gravitational first order phase transition". It has to be noted, however, that contrary to the 
liquid/vapor transition, the two phases cannot coexist in the present situation. 

For A t (fi) < A < A c , the non degenerate solutions are metastable but we may suspect 
|| that they are long-lived so that they are physical. These solutions are insensitive to the 
small-scale regularization and depend only on the long-range gravitational interaction. In the 
limit /i — > +oo, the transition energy A 4 (/x) goes to — oo and we recover the classical spiral 
of Fig. [I]. This spiral is formed by the metastable states of the "gaseous" phase (points A). 
The "condensed" phase (points C) is superposed to the A-axis. These states have an infinite 
central density and an infinite temperature. The unstable branch (points B) coincides with the 
spiral but these states physically differ from the "gaseous" states (points A) by the presence 
of an infinitesimal "germ" with negligible mass and energy. Therefore, what we see actually 
in the limit fi — > +oo are the metastable states. It is plausible that these metastable states 
will be selected by the dynamics (on relevant time scales) even if states with higher entropy 
exist. In fact, depending on its topology (i.e., the form of the profile) an initial condition with 
A > A t (n) can either relax towards the local entropy maximum (gaseous phase) or collapse 
towards the global entropy maximum (condensed phase). Alternatively, for A < A t (fi) an 
initially "condensed" configuration can remain frozen in this metastable state or explose in a 
"gaseous" state with more entropy. Therefore, the choice of a stable equilibrium state does not 
only depend whether the equilibrium solution is a local or a global entropy maximum. It is 
more dependant whether the initial condition lies in the "basin of attraction" of the equilibrium 
state or not. The characterization of this basin of attraction requires a nonequilibrium analysis 
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Figure 7: Equilibrium phase diagram for Fermi-Dirac spheres with a degeneracy parameter 
/x = 10 5 . Points A form the "gaseous" phase. They are global entropy maxima (GEM) for 
A < A t (/i) and local entropy maxima (LEM) for A > A t (/x). Points C form the "condensed" 
phase. They are LEM for A < A t (/x) and GEM for A > A t (/x). Points B are unstable saddle 
points (SP) and contain a "germ" . This figure exhibits in particular a first order phase transition 
in the microcanonical ensemble. 




Figure 8: Density contrast as a function of energy for self-gravitating fermions (/x = 10 5 ). This 
figure can be compared with Fig. |2] for classical isothermal spheres. Points A (gaseous phase) 
have a low density contrast. Points B and C contain a central nucleus with high density |J. 
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Figure 9: Entropy of each phase versus energy for fi = 10 5 . A phase transition occurs at Af(/x) 
at which the two stable branches (solutions A and C) intersect. The unstable solutions B always 
have smaller entropy. 

which is not attempted in the present paper. A first step in that direction was performed by 



Youngkins & Miller [jllj by using a one- dimensional spherical shell model and by Chavanis et al. 
pBj with the aid of a simple relaxation equation derived from a maximum entropy production 
principle [pl|. These preliminary works reveal that the structure of this basin of attraction is 
extremely complex so that the final state of the system cannot be easily predicted from the 
initial condition when several equilibrium states exist. In addition, the structure of this basin of 
attraction depends whether the system is described by the microcanonical or by the canonical 
ensemble 0, 0. 

However, for A > A c the metastable phase completely disapears and, in that case, the system 
must necessarily collapse. This transition is associated to what has been traditionally called 
the "gravothermal catastrophe" || in the case of classical point masses. For systems described 
by the Fermi-Dirac statistics, the core ultimately ceases to shrink when it becomes degenerate. 
In that case, the system falls on to the global entropy maximum (points D) which is the true 
equilibrium state for these systems. This global entropy maximum has a "core-halo" structure, 
with a degenerate core and a non degenerate halo. This phase transition is sometimes called a 
zeroth order phase transition [^] since it is associated with a discontinuous jump of entropy 
(in the classical limit, the entropy of the condensed phase is infinite). In fact, this does not 
correspond to a true phase transition (in the usual sense), not even to an instability but simply 
to the sudden disapearance of the "gaseous" phase. It has to be noted that the degenerate 
nucleus resulting from this gravitational collapse has a relatively important mass and a very 
small radius (for /i = 10 5 and A ~ A c , we have typically M*/M ~ 0.22 and R*/R ~ 5 1CT 3 ). 
This massive nucleus ("fermion star") can have important astrophysical implications and, in 
the context of dark matter, may mimick the effect of a black hole at the center of galaxies (see 
Sec. 



12 



3.3 Small values of the degeneracy parameter (/i = 10 3 ) 

When the degeneracy parameter is sufficiently small, there exists only one critical point of 
entropy for each value of energy (see Fig. [10]) and it is a global entropy maximum. Therefore, a 
sufficiently strong degeneracy suppresses the phase transitions in the microcanonical ensemble, 
including the "gravothermal catastrophe" (see Fig. [TT| ). For high energies (small values of A) 
the solutions almost coincide with the classical isothermal spheres. When the energy is lowered 
(large values of A) the solutions take a "core-halo" structure with a partially degenerate nucleus 
surrounded by a dilute Maxwellian atmosphere. It is now possible to overcome the critical 
energy A c = 0.335 and the critical density contrast TZ = 709 found by Antonov || for classical 
particles. In that region, the specific heat is negative, which is allowed in the microcanonical 
ensemble. As energy decreases further, more and more mass is concentrated into the nucleus 
(which becomes more and more degenerate) until a minimum accessible energy, corresponding 
to A max (/j,), at which the nucleus contains all the mass. In that case, the atmosphere has been 
"swallowed" and the system has the same structure as a cold white dwarf star |14}| . This is a 



relatively singular limit since the density drops to zero at a finite radius whereas for partially 
degenerate systems, the density decays like r~ 2 at large distances. We can study the formation 
of this compact object by defining an order parameter k = R95/R, where -R95 is the radius 
of the sphere which contains 95% of the mass fl0| . This parameter is ploted as a function of 



energy in Fig. and the diagram is similar to the one obtained by Follana & Laliena fll0| 
with a different regularization of the potential at the origin. For high energies, the density 
varies smoothly with the distance and k ~ 1. For low energies, the system spontaneously forms 
a dense core containing more and more mass so that k ^ k c « 1 (with k c ~ 6.678/yU 2 / 3 , 
estimated from the mass-radius relation ( |3"3] ) of a completely degenerate nucleus). We observe 
that the order parameter varies rapidly in the region of negative specific heats but remains 
continuous. According to Cerruti-Sola et al. P5[| , this is the mark of a second order phase 
transition at A = A (corresponding to the point of minimum temperature r] c ). At that point, 
the specific heat is infinite and turns from positive to negative; more precisely, C = dE/dT 
diverges like (A — Ao) _1 and ±(r/ c — r])~ l l 2 [25]. In fact, this "second order phase transition" is 
not really a phase transition; it just corresponds to the "clustering" of the self-gravitating gas 
as its energy is progressively reduced. 

If we now consider the canonical situation, we see that several solutions exist at the same 
temperature. A first order phase transition occurs at T)t(^) an separates a "gaseous" phase 
from a "condensed" phase. The interpretation is the same as in the microcanonical ensemble. 
The energy is discontinuous at the transition so that a large amount of latent heat is released. 
The specific heat is also discontinuous but remains positive. For 77 > the "gaseous" states 
are metastable but they probably are physical. This metastable branch completely disapears 
at T) = 7] c and the system undergoes an "isothermal collapse". This phase transition is more 
radical than the previous one since it is marked by the disapearance of the metastable phase. 
For self-gravitating fermions the "isothermal collapse" ends up on a compact state with a core- 
halo structure (for \i = 10 3 , M* ~ M and R* ~ 6.7 10~ 2 i? at 77 ~ rj c ). Since the free energy is 
discontinuous (see Fig. O), this could be called a zeroth order phase transition PS[. It has been 



shown in Ref. |25j that the point of minimum temperature 7] c coincides with the Jeans instability 



criterion. More precisely, the condition of instability of an isothermal gas sphere with respect to 
linear perturbations of the Navier-Stokes equation can be written R > (r] c /3) 1 ^ 2 Lj ~ Lj, where 
Lj is the Jeans length. Since r] t < i] c , a phase transition (corresponding to a nonlinear evolution 
of the system) can occur at scales much smaller than the Jeans scale (R = (rjt/SY^Lj <C Lj). 
This might explain the formation of smaller objects than usually achieved with the ordinary 
Jeans instability. This idea has been developed by Stahl et al. M in relation with planet 
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Figure 10: Equilibrium phase diagram for Fermi-Dirac spheres with a degeneracy parameter 
/j, = 10 3 . Points A form the "gaseous" phase. They are global maxima of free energy (GFEM) 
for f] < i] t and local maxima of free energy (LFEM) for r\> r\t- The reverse is true for points 
C in the "condensed" phase. Points B are unstable saddle points (SP). This figure exhibits 
in particular a first order phase transition in the canonical ensemble. The equality of the free 
energy of the two phases (J a = Jc) at the transition temperature r\ t implies the equality of 
the areas delimited by the curve and the plateau, i.e. J A (v — r) t )dA = 0. This is similar to 
Maxwell's construction in the theory of the van der Waals gas. 
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Figure 11: Equilibrium phase diagram for self-gravitating fermions with different values of the 
degeneracy parameter //. The zeroth and first order phase transitions are suppressed in the 
microcanonical ensemble for /i < 2600 and in the canonical ensemble for /i < 82.5. For large 
values of /i, the curve makes several rotations before unwinding. The criterion of Katz tells 
us that one mode of stability is lost each time the curve rotates clockwise and regained as the 
curve rotates anticlockwise. Therefore, only the upper and lower branches are entropy maxima. 
For fi — > +oo (classical limit), the curve winds up indefinitely and tends to the spiral of Fig. [1] 
as discussed in the text. 




Figure 12: Evolution of the order parameter with the energy for \x = 10 3 . This figure illustrates 
the "clustering" of the self-gravitating gas as energy is lowered. The presence of the "bump" 
at A ~ A c was previously noted by Follana & Laliena flM in their model. 
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Figure 13: Free energy of each phase versus inverse temperature for /i = 10 3 . The interpreata- 
tion is the same as in Fig. |9|. 

formation. However, this picture is limited by the realization that the evolution of the system 
depends on a complicated notion of "basin of attraction" which is difficult to describe in detail 
(see Sec. |3.2|) . In particular, it is not clear whether the true phase transition will occur at r\ t . 



4 Tricritical points 

We have indicated in the preceding section that the phase transition in the microcanonical en- 
semble disapears at a critical degeneracy parameter Hmtp — 2600 and that the phase transition 
in the canonical ensemble disapears for fi > [Ictp — 82.5. These points at which the A — i] 
curve presents an inflexion are called tricritical points (TP) in the langage of phase transitions. 



In Fig. [14|, we have enlarged the phase diagram near the tricritical point in the canoni- 
cal ensemble. It is located at Aqtp — 0.5 and tjctp — 3.06. In the canonical ensemble, the 
oscillations of the A — rj curve for /x > ^ctp are replaced by a horizontal Maxwell plateau con- 
necting the gaseous phase (left) to the condensed phase (right). This characterizes a canonical 
first order phase transition at f]t{^)- This diagram exhibits a close analogy with the classical 
gas/liquid transition, the liquid phase beeing the counterpart of the (gravitational) Fermi con- 
densate. At the tricritical point CTP, the plateau disapears and the specific heats diverges like 
C ~ (A — A c )~ 2 ~ \r/ — i] c \~ 2 / 3 . This is similar to the A-transition in 4 H e . Therefore, the first 
order phase transition becomes second order at the canonical tricritical point. 

The gravitational Fermi gas diagram is nevertheless more complex than the liquid/gas dia- 
gram because it presents another tricritical point MTP in the microcanonical ensemble. In Fig. 



15] , we have enlarged the phase diagram near this tricritical point. It is located at Amtp — 0.38, 
rjMTP — 1-68. The interpretation is the same as in the canonical ensemble except that the 
plateau is now vertical as it corresponds to a discontinuity of temperature at a transition en- 
ergy A 4 (/i) (microcanonical first order phase transition). We have denoted by r} ga s and rj con d the 
values of the inverse temperature of the two phases at the transition energy. At the tricritical 
point r] gas = r] cond . We have also indicated in the figure the point of minimum temperature r\ c 
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Figure 14: Enlargement of the phase diagram near the tricritical point in the canonical ensem- 
ble. The Maxwell construction determining the transition temperature r) t (ii) is done explicitly 
(dashed areas). For [ictp = 82.5 the Maxwell plateau disapears and the A — rj curve presents 
an inflexion point at Actp — 0.5, tjctp — 3.06. At that point, the specific heat becomes infi- 
nite and the transition is second order. This diagram is remarkably similar to the liquid/gas 
transition for an ordinary fluid. 




Figure 15: Enlargement of the phase diagram near the tricritical point in the microcanonical 
ensemble. The Maxwell construction determining the transition energy A t (//) is done explicitly 
(dashed areas). For /imtp = 2600 the Maxwell plateau disapears and the A — rj curve presents 
an inflexion point at A M tp — 0.38, t]mtp — 1-68. We have indicated different characteristic 
temperatures, as described in the text. 
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at which the specific heat diverges (second order microcanonical phase transition). 

4 i 1 1 1 1 1 1 1 




Figure 16: Phase diagram of the self-gravitating Fermi gas in the \i — r\ plane. The solid line 
gives the transition temperature at the canonical first order phase transition. The dashed lines 
give the temperatures of the two phases (gas and condensate) at the microcanonical first order 
phase transition. The dotted line gives the critical temperature 7] c at the microcanonical second 
order phase transition. CTP and MTP are the tricritical points in the two ensembles. 



In Fig. 16, we have ploted the values of the transition temperature rj t in the canonical 
ensemble and the values of the characteristic temperatures rj gas and r) con d in the microcanonical 
ensemble as a function of the degeneracy parameter \i. These curves characterize first order 
phase transitions in the canonical and microcanonical ensembles. We have also represented the 
minimum temperature r] c marking the onset of a microcanonical second order phase transition. 
At the canonical tricritical point CTP, this branch connects the ^-branch. This is consistent 
with the appearance of an isolated second order phase transition in the canonical ensemble 



at CTP. This diagram shares some similitudes with the one obtained by Barre et al. [12] in 
their analysis of the Blume-Emery-Griffiths (BEG) model with infinite range interactions. This 
system also displays an inequivalence of ensembles, regions of negative specific heats and two 
tricritical points (one in each ensemble). This analogy suggests that these properties are related 
to the long-range nature of the interactions more than to the details of the model. This implies 
a kind of universality for this type of systems. However, there are also noticable differences 
between the two models. In particular, in the BEG model, second order phase transitions 
are characterized by a discontinuity of the specific heat (angular point) while in our model of 
self-gravitating fermions the specific heat is always continuous except at the critical point r] c 
at which it diverges (this is our definition of a second order phase transition). Therefore, there 
still exists a second order phase transition above the canonical tricritical point in the BEG 
model but not in the self-gravitating Fermi gas. Another consequence of the continuity of the 
specific heats in our model is that the second order critical line r] c — /i (dotted line) does not 
connect the first order critical lines (dashed lines) at the microcanonical tricritical point MTP 
but slightly after, unlike in the BEG model fL2|. 
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5 A simple analytical model for self-gravitating fermions 



The previous study has revealed that self-gravitating fermions can undergo a phase transition 
from a "gaseous" phase to a "condensed" phase. We shall now propose a simple analytical 
model to describe this phase transition more conveniently. As we shall see, our model can 
reproduce remarkably well the essential features of the equilibrium phase diagram and it can 
be used to determine the dependance of the critical parameters with the degeneracy parameter 



5.1 The "gaseous" phase 

The "gaseous" phase can be represented by a homogeneous distribution of particles with a 
Maxwellian distribution of velocities 

(24) ^WrfT^- 

The relation between the energy and the temperature is given by 

(25) E = hlT - 
and the entropy by 

VoS 3 , . mN / 3M \ , 3 

26 — = - In 2ttT - n + n 77 + -. 

K J M 2 \AnR 3 J /0 2 



Introducing the normalized energy A and the normalized temperature rj defined in section 2J2 
we can rewrite the previous equations in the form 

(27) A = - - — , 

K ' 5 277' 



(28) 



- = — - In 77 + In a + — H — In tt — In 6. 
M 2 ' p 2 2 



These equations correctly describe the gaseous phase for high energies and high temperatures 
(i.e. low density contrasts). Of course, it can not reproduce the spiral behaviour of the classical 
phase diagram which is an intrinsic property of the Emden equation. An analytical expression 
for this spiral has been given in Ref. 



25|j in the asymptotic limit of high density contrasts. 



5.2 The "condensed" phase 

For the "condensed" phase, we shall improve the core-halo model proposed by Chavanis & 
Sommeria 0. We assume that the "core" is completely degenerate and we denote by M*, R* 
and its mass, radius and energy. In that limit, the distribution function is a step function: 
/ = 7/0 for v < v max and / = for v > v max . Of course, for a self-gravitating system, the 
maximum velocity v max is a function of the position. For this simple distribution function, the 
pressure and the density are given by 

(29) p J/„,Vv = ^% 
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(30) p = J fd\ = 4^^. 

Eliminating the velocity between these two relations, we find that the equation of state of a 
completely degenerate system is that of a polytrope with index 7 = 5/3 (or n — 3/2), 

< 31 > K 4 (4) V3 

These results are of course well-known from the theory of white dwarf stars [13] and they are 
repeated only in order to determine the constant K in the present context. Now, from the 
theory of polytropic spheres |14|], the mass-radius relation is given in the general case by 

(32) K = N n GM {n - 1)/n R {3 - n)/n , 

where N n is a constant depending on the index of the polytrope. For n = 3/2, one has 
N3/2 = 0.42422.... Therefore, the relation between the mass and the radius of our degenerate 
nucleus is 

(33) M*Rl = with X - 5.9723 10" 3 . 

r] G 



On the other hand, its energy is given by 

3 GMl 



(34) E, 



7 iL 



and its entropy is equal to zero since the distribution function is unmixed (/ = 770). 

By shrinking, the nucleus releases an enormous amount of energy which heats the envelope. 
The envelope behaves therefore like an ordinary gas maintained by the walls of the box so that 
its density is approximately uniform. Its energy and entropy are therefore given by 

(35) S 3(M-M,)T 3GM,(M-M.) 3G(M - M, ? 



2 y ' 2R 5R 



(36) 



Vo S = (M — Af„) 



3 3 
- ln(27rT) - ln(M - M*) + In V + In r] + - 



Contrary to our previous paper ||, we have not neglected the potential energy of the envelope 
as compared to its thermal energy. This sensibly improves the agreement with the full numerical 
solution. However, we have still assumed that the core is much smaller than the halo, so that 
V = |vri? 3 represents the total volume of the system. For calculating the potential energy, we 
have used the formula 



W = -4ttG / pM(r)rdr 



(37) 



where M(r) is the mass contained within the sphere of radius r. This formula is valid for an 
arbitrary spherically symmetrical distribution of matter |22 . 

Adding Eqs. fl3"4"D (|3li| ) and expressing the radius of the core as a function of its mass, using 
Eq. (|33D, the total energy of the system is given by 

(38 ) E = + \(M - M,)T - MM,(M-M,) _ ^M-M^ 
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For a given value of energy (microcanonical ensemble), this relation determines the temperature 
T as a function of the core mass. Therefore, the entropy ( |3"6| ) is a function of M* alone. The 
mass of the nucleus at equilibrium is determined by maximizing the entropy with respect to 
M*. After simplification, the condition dS/dM* = is found to be equivalent to 



ln(M - M») 



3G 
2RT 



(M - 2M*) 



6G 
5RT 



{M-M* 



-InT 

2 



X 



1/3 T 



(39) 



In r]o + In V - 1 + - ln(27r) 



We obtain the same relation by maximizing the free energy J at fixed mass and temperature. 
Eqs. (|38| ) - (|39|) completely determine the equilibrium phase diagram of self-gravitating fermions 
in the framework of our analytical model. For convenience, we shall re-express these equations 
in a dimensionless form. To that purpose, we introduce the fraction of mass a contained in the 
core such that M* = ocM. In terms of a, the radius of the core is R*/R = 6.678 a^ 1 ^ 3 ^ 2 ^ 3 
(see Eq. (]3~3f)). On the other hand, we set 



(40) 



A 



1 



0.149756..., 



C 



1 



In7r-ln6 - 1 = -2.21939. 



Introducing furthermore the dimensionless energy A and the dimensionless temperature r\ de- 
fined in Sec. the equations of the problem become 



(41) 



A=^A/i 2 /V/ 3 --(l 
7 P 2 V 



a) 



-a(l — a) 



'I -a) 



(42) 



ln(l - a) + - In 77 + A/i 2/3 



a 



4/3 



— a 



In fx + C, 



(43) 



VqS 
M 



'I -a) 



— In rj — ln(l — a) + In fi + C 



(44) 



M 



M 



+ r]A. 



The equilibrium phase diagram is represented in Fig. |17|for a degeneracy parameter fi = 10 5 . 
It provides a fairly good agreement with the full numerical solution of Sec. \i.2\ (see Fig. 
0). Of course, we cannot expect to reproduce exactly the numerical results in view of the 
approximations made in our analytical model. In particular, except for very low energies, 
the core is only partially degenerate and this is responsible for the quantitative discrepencies 
observed between the two diagrams. However, the qualitative behaviour is the same and this 
is essentially what was attempted by our analytical approach. In particular, we recover the 
three types of solutions previously studied. The solutions on the upper branch (points A) form 
the "gaseous" phase. They can also be considered as a particular limit of the core-halo model 
with a = 0. The solutions on the lower branch (points C) form the "condensed" phase and 
the solutions on the intermediate branch (points B) are similar to the "gaseous" states (points 
A) except that they possess a small central nucleus (a "germ"). To determine the stability of 



these solutions we have ploted in Fig. [L8| their entropy as a function of the core mass a for 
a given energy. The curve S(a) has the usual 'W-shape" characteristic of phase transitions. 
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Figure 17: Equilibrium phase diagram obtained from our analytical model (/i = 10 5 ). It 
compares relatively well with the full numerical solution reported in Fig. 




Figure 18: Entropy as a function of the mass of the nucleus for /i = 10 5 and A = —0.1 > A t . This 
curve is obtained by using the analytical formulae (f4T|) ([42]) ([43|) . For A > A t , the "condensed" 
state (point C) is a global entropy maximum and the "gaseous" state (point A) a local entropy 
maximum. The solution with the "germ" (point B) is an entropy minimum. The fraction of 
mass contained in the nucleus is relatively large for /z = 10 5 (see Fig. |T9|) but it decreases as 
the classical limit is approached. As /i — > +oo, a con d — > 0, a germ — > with a con d 3> a germ . 
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Figure 19: Evolution of the fraction of mass contained in the nucleus as a function of energy 
for a degeneracy parameter \l = 10 5 . 
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Figure 20: Equilibrium phase diagram obtained from our analytical model (/i = 10 3 ). It 
compares relatively well with the full numerical solution reported in Fig. [TO. 



23 




Figure 21: Free energy as a function of the mass of the nucleus for \x = 10 3 and rj = 0.74 < r] t . 
For T] < T]t, the "gaseous" state is a global maximum of free energy and the "condensed" state 
is a local maximum of free energy The solution with the "germ" is a minimum of free energy. 
The fraction of mass contained in the condensate is close to one (see Fig. |22[). In the limit 
jj, — > oo, otcond — > 1 and a germ — > 0. 




Figure 22: Evolution of the fraction of mass contained in the nucleus as a function of the inverse 
temperature for a degeneracy parameter \i = 10 3 . 
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The "gaseous" state (a = 0) can be considered as an entropy maximum although it does not 
correspond to the condition dS/da = 0. The other entropy maximum corresponds to the 
"condensed" state (point C) and the entropy minimum to the solution with the germ (point 
B). For A > A t , the "condensed" states C are global entropy maxima and the "gaseous" states 
A are local entropy maxima (the reverse is true for A < A t ). However, to pass from A to C, we 
have to cross an entropic barrier constituted by the solution B. This requires that the entropy 
decreases which is not possible for an isolated system. Therefore, depending whether the initial 
fraction of mass ao is smaller or larger than a germ , the system can either relax towards the 
"gaseous" state or collapse towards the "condensed" state. Of course, this argument assumes 
that a is the only degree of freedom in the system, which is clearly an idealization. As mentioned 
previously, the real "basin of attraction" is much more complicated. 

The equilibrium phase diagram corresponding to a degeneracy parameter /i = 10 3 is rep- 
resented in Fig. |20] and it compares relatively well with the full numerical solution (see Fig. 



10). The minimum energy, corresponding to A max (/j,) } is reached when all the mass is in the 
degenerate core at zero temperature (a — > 1, T] — > +oo). Eq. (|4l|) then yields 

(45) A ma M = ^A/i 2 / 3 . 
For A — > A max (fi) and for sufficiently large n, one has 

(46) 1 - a ~ e- A ^\ A max - A ~ X^e~ x ^. 

This shows that the atmosphere is swallowed exponentially rapidly when we approach the 
minimum energy. Fig. [2t] displays a clear phase transition in the canonical ensemble. The 
value of the temperature of transition is very close to the exact value found in the numerical 
approach. In the present case, we can safely consider that the core of the condensate is com- 
pletely degenerate so that the quantitative agreement with the exact solution is better than 
in the microcanonical ensemble. Once again, the stability of the equilibrium states can be 
determined by considering the variation of the free energy with the core mass (see Fig. |2T|) . 
For 7] < rj t (resp. rj > i] t ), the "gaseous" state is a global (resp. local) maximum of free energy 
and the "condensed" state a local (resp. global) one. The description of the phase transition 
is the same as in the microcanonical ensemble. 



5.3 Scaling laws in the limit ji 



-oc 



We shall now determine the behaviour of the critical parameters as [l — > +oo (classical limit). 
In the microcanonical ensemble, the phase transition occurs close to the maximum energy E* of 
the condensed phase, corresponding to A*(/i). Using Eqs. (|4]) fl42|) , we find that the condition 
of energy maximum dA = is equivalent to 

2 o 



(47) 



+ A/i 2//3 



a 



4/3 



9 

TJ + -71 



a 



a) 



-X^a^r] - - J— - -77 
3 1 — a 5 



This equation, together with Eqs. ((y]) (|42|) , determines the function A*(/x). Now, in the limit 
\l — ► +00, the fraction of mass contained in the nucleus goes to zero while the temperature 
increases. Taking the limit /i — > +00, a —* and 7] — > in Eq. (|47|) , we obtain the relation 

(48) Xfi 2 / 3 a 7 / 3 V = 2. 

Considering the limiting form of Eqs. (fll]) ( |4~2| ) in the same approximation, we find 

-A/^V/ 3 - — , 
7 P 277' 



(49) 



A 
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(50) | InTy + A/i 2/3 a 4/3 r/ = In /x. 

Therefore, the thermodynamical parameters at the point of maximum energy behave with the 
degeneracy parameter like 

1 (ln/i) 7 / 3 ^2/3 

(51) a* ~ , r/* ~ 7; , A, - 



ln/x' * /i 2 / 3 ' (ln/i) 7 / 3 

The transition point A t (fi) is determined by equating the entropy of the two phases. This 
yields 



(52) (l-«) 



3 5 
• ^ ln ^cond - ln(l - a) + In fi + C + - 



3 5 
--\nr] gas + ln/i + C + -, 



where is the temperature of the condensed phase, rj gas the temperature of the gaseous 
phase and a the mass contained in the condensate at the transition point. Considering the 
limiting form of Eqs. fl27|) ( f4~T|) (^) when /i — ■> +oo, a — > and ?7 — > 0, we find that 

(53) A = — = h^a 1 ' 3 - — ^— 



^Vgas 7 ^Vcond 

3 

(54) - ln r] cond + Xfi 2/3 a i/3 r] cond = ln /x. 
Solving for t] cond and 7/ 5as , we get 

2 ln/i 1 A/i 2 / 3 a 4 / 3 / 4 , , 

(55) ricond ~ — l--aln/i . 



A/i 2/3 a 4/ 3 ' %m 2 ln/i V 7 



Substituting these results in Eq. (p2|) , we find that the mass contained in the condensate 
behaves like 

(56) a 



ln /i 

Combining the foregoing relations, we find that the energy of transition is given by 

u 2/3 

(57) A, ~ 



(ln/i) 7 / 3 ' 

In the classical limit, At — ► — oo, so that the "gaseous" states are always metastable as previ- 
ously discussed. The temperature of the two phases at the transition point behave like 

(ln/i) 7 / 3 

(58) r] gas ~ r] cond ^ 2/3 ■ 

The jump of temperature is given by 

1 1 /i 2 / 3 



(59) 



Vcond Vgas (ln/i) 7 / 3 ' 



Finally, we can study the behaviour of the fraction of mass contained in the germ and in the 
condensate as \i — > +oo for a given value of energy A. For the condensate (points C), a — > 
and r) — > so that Eqs. ([H]) fl4"2] ) simplify in 

(60) A = ^/v/3, 
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(61) ^lnr/ + A y u 2/3 a 4/3 r/ = ln/i. 

The first equality simply means that the potential energy of the core tends to — oo so that the 
temperature of the halo must rise to +00 so as to maintain the total energy fixed. Combining 
Eqs. (|60D (l5l]) and using Eq. (f43f ), we obtain 



1 (IM 7 / 3 q , 

(62) Oicond ~ -, , Vcond ~ ^ , ^cond ~ In /I. 

In fj, \i l l 6 

Since the entropy diverges as fi — > +00, we recover the well-known fact that there is no global 
entropy maximum for a classical self-gravitating gas. Note that the divergence of entropy is 
relatively slow (logarithmic). For the germ (points B), we have in the classical limit 



(63) Otgerm * 0, T]germ * ^Igasy Sg erm ► S, 



gas- 



These results simply reflect the fact that the unstable branch approaches the gaseous branch 
as n — > +00 but still differs from it by a presence of a small germ. The size of the germ is 
determined by Eq. (|42|). For /i — > +00, it leads to 



(64) X^ 3 a^\ as = \np 

so that 

(ln^u) 3 / 4 



(65) 



germ 1/2 • 

For a given energy, the mass contained in the condensate and in the germ goes to zero as we 
approach the classical limit. However, the relation 



Otgerm (h' AO 



7/4 



(66) ^HL ~ « 1, 

®-cond 

indicates that the size of the germ is much smaller than the size of the condensate. 



These results are consistent with the proof given in Ref. for the absence of global 
entropy maximum in the microcanonical ensemble. In words, we can make the entropy diverge 
by approaching an arbitrarily small fraction of particles in the core (M# <C M) so that the 
potential energy goes to —00. Since the total energy is conserved, the temperature must rise 
to +00 and this leads to a logarithmic divergence of the entropy. This "natural" evolution is 
confirmed by dynamical models of self-gravitating systems (see the discussion in Ref. [p6| ). It 



is found that the gravitational collapse of classical point masses leads to a finite time singularity 
(the central density becomes infinite in a finite time t co u) with a slow algebraic divergence of 
the temperature and a logarithmic divergence of the entropy as t — > t co u. In addition, the mass 
contained in the core tends to zero as t — > t co u (the density profile at t — t co u is close to a power 
law p ~ r~ a , with a ~ 2.2 < 3) in agreement with our previous observations. 

We now consider the canonical situation. In the limit \i — > +00, the phase transition occurs 
close to maximum temperature of the condensed phase, corresponding to 77* (/V). Using Eq. 
fl42|), we find that the condition of temperature maximum drj = is equivalent to 



(67) % 2/3 a^ V -^-- 9 V = 0. 

3 1 — a 5 
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This equation, together with Eq. fl42|), determines the function 77* (/i.) . In the limit /i — > +00, the 
fraction of mass contained in the condensate comes close to unity and the maximum temperature 
increases. Taking the limit a — > 1 and rj — > in Eq. fl6"7j), we get 



(68) 



4 



1 — a 



Simplifying Eqs. (|4l|) (|42|) in the same approximation, we obtain 



(69) 



A = -A/x 2 / 3 
7 



5(i 

2^ 



a)-, 
V 



(70) 



ln(l — a) + - In r\ + A/i 2//3 7/ = In \i. 



We find therefore that the parameters at the point of maximum temperature behave like 



(71) 



In /1 



ln/x' 



2/3' 



A* A maa ,(/z) 



2/3 



(hi/i) : 



The transition point rjt (/i) in the canonical ensemble is obtained by equating the free energy 
of the two phases. Using the results of Sec. |5.1| and Sec. |5.2|, we obtain the general relation 



3 3 

- In 7/ + ln/i + C + 1 + -77 = (1 - a) 



3 5 
-- In T] - ln(l -a)+\nfi + C + - 



(72) 



+-\r]^ 3 a 7/3 - -(1 - a) + -770(1 - a) + -7/(1 - «) 2 



where a is the mass contained in the condensate at the transition point. In the limit /j — > +00 
and a — > 1, this relation reduces to 



(73) 
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The free energy of the gaseous phase (l.h.s.) is dominated by the contribution of the entropy 
and the free energy of the condensed phase (r.h.s.) is dominated by the energy of the core. We 
find therefore that the temperature of transition behaves like 



(74) 
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The mass contained in the condensate is determined by Eq. (f|2). In the limit /i — > +00, it 
simplifies in 



(75) 



ln(l — a) + — In rj + A/i 2 ^ 3 7/ = In /1. 



Solving for a in Eqs. (|T3"D (fT5|), we get 

(76) 1 - a 
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According to Eqs. (|27|)([4lD ([74]) , the energy of the gaseous phase and the energy of the con- 
densed phase behave at the transition point like 



(77) 
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The jump of energy (latent heat) is 



(78) A cond - A gas ~ fi 2/3 ( 1 + 



4m/^ 

Finally, let us determine the fraction of mass contained in the condensate and in the germ 
as fi — > +00 for a given value of temperature. For the condensate (points C), a — > 1 so that 
Eq. (g2p leads to 

(79) l-a cond ~e- x ™ 2/ \ 
The free energy of the condensed state behaves like 

(80) J ~ KnaxilAV ~ 

The divergence of the free energy in the canonical ensemble is more rapid than the divergence of 
the entropy in the microcanonical ensemble. This is simply because the free energy is dominated 
by the divergence of the (potential) energy while the entropy is dominated by the divergence 
of the logarithm of the temperature. For the germ (points B), we have in the classical limit 

(81) CVgerm * A > Ag as J ► Jg as . 



Using Eq. (§2]), we find that 

(82) a. - <'^ 3/1 



<-germ 1/2 • 

In the limit fi — ► +00, the germ contains almost no mass while the condensate contains almost 
all the mass. 

These results again agree with the proof given in Ref. [EH] for the absence of a global max- 
imum of free energy in the canonical ensemble. In words, we can make the free energy diverge 
by collapsing the mass M to a point (the divergence of J is simply due to the divergence of the 
potential energy; the entropy has a weak, logarithmic, negative divergence). This argument 
would suggest that the natural evolution of a system of classical point masses in the canonical 
ensemble is to develop a density profile in the form of a ^-function with all the mass at r = 0. 
This is not what numerical simulations of gravitational collapse show (see the discussion in Ref. 
2T)]). It is usually found that when the system is held at a fixed temperature, the self-similar 



collapse leads to a density profile close to the power law p ~ r~ 2 at t — t co u. This profile has a 
vanishing mass at r = and its free energy is not divergent. Therefore, a finite time singularity 
prevents the system from reaching arbitrarily large values of the free energy . It is not known 
whether other solutions of these dynamical equations (not necessarily self-similar) can lead to 
the expected 5-function with J = +00. 



6 Classical gas with a short distance cut-off 

The case of a self-gravitating gas with a short distance cut-off was first considered by Aronson 
& Hansen || and more recently by Stahl et al. (see discussion in Sec. |6^) . The equilibrium 



phase diagram of this system is similar to the one obtained for self-gravitating fermions (see 
Figs. |7]-|T0"1). Indeed, the degeneracy parameter \i plays the same role as the inverse of the 
short distance cut-off a. The interpretation of the phase transitions is therefore similar but 
the dependance of the critical parameters with the cut-off is different. We shall therefore 
reformulate our analytical model to the case of a classical hard spheres gas and determine how 
the previous results are modified in this new situation. 



29 



6.1 The "gaseous" phase 



We model the gaseous phase by a uniform distribution of matter occupying the whole container. 
The energy and the entropy are therefore given by 



(83) 
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In dimensionless variables, these equations can be rewritten 
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~ 5 2r/' 



In the expression ([ 
(87) 

which plays no role in the problem. 



S/N = --\nr]. 
for the entropy, we have not written the constant term 
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6.2 The "condensed" phase 

We model the "condensed" phase by a nucleus and an atmosphere, each of uniform density. 
The velocity distribution of the particles is assumed to be Maxwellian with temperature T. Let 
R* be the radius of the nucleus and Af* the number of particles that it contains. We introduce 
an excluded volume ~ ^na 3 around each particle, where a can be regarded as the "effective" 
size of the particles. R* and A^* are therefore related to each other by a relation of the form 

Rl = 4gN*a\ 



where g is a geometrical factor with order of magnitude unity which depends on the nature of 
the close-packing (see Ref. ||). The energy of the core is 



(89) 

and its entropy 
(90) 

For the halo, we have 
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Like in the case of fermions, we have considered that the volume of the nucleus is much smaller 
than the volume of the halo. Adding these expressions and using Eq. fl8"8|) to express the radius 
of the core as a function of its mass M* = N*m, we obtain for the whole configuration 



(93) 
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Let a = M*/M denote the fraction of mass contained in the nucleus. We also introduce the 
filling factor 



(95) 



R 



aiigN) 1 / 3 ' 



which can be regarded as an inverse normalized hard spheres radius. The case of point masses 
corresponds to the limit « — > +oo. Clearly, « plays the same role as the degeneracy parameter 
in section E3| In dimensionless form, the equations of the problem become 



(96) 



A = -— + \iia^ + -a(l -a) + -(1 
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(97) S/N = -- In T) - 3a lnu - (1 - a) ln(l - a), 

where we have again eliminated the constant ([8j]) from the expression of the entropy. We now 
determine the mass of the nucleus by maximizing the entropy ( |9~TD at fixed energy. This yields 
the relation 

(98) 1 + ln(l - a) - 3 In fi = + ^ar] - fir]a 2/3 . 

Eqs. (pop (p7|) (|9~8|) determine the equilibrium phase diagram of a classical hard spheres gas in 
the framework of our analytical model. The description of this diagram (and its dependance 
on the parameter /x) is similar to the one given in section |5] for self-gravitating fermions. The 
minimum energy corresponds to the configuration for which all the mass is in the nucleus at 
zero temperature. Taking the limit a — > 1 and 77 — > +00 in Eq. (p6|), we get 

3 

(99) Knax{lj) = TV- 

5 

For A — > A max and for sufficiently large values of /z, Eqs. (|96l) (|98l) yield 

(100) 1 - a ~ e- w , A max — A~^-. 

Note that the relation between the temperature and the energy is different from the corre- 
sponding one for self-gravitating fermions (see Eq. (|46|) ). 
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6.3 Scaling laws in the limit /i —> +00 

The derivation of the scaling laws for the critical parameters of a classical hard spheres gas is 
essentially the same as for self-gravitating fermions (section |5.3| ). We shall directly give the 
results without detailed discussion. 

In the microcanonical ensemble, the point of maximum energy, corresponding to A*(//), is 
determined by the relation 

/ x 2 9 / 3 9 9/ o\ 2 2 1/, 9 1 

(101) —77 I — — —a + jJLOt ' ] = -^a 1 ~ 
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In the limit /t — > +00, the thermodyanical parameters behave at the point of maximum energy 
like 

/im , 1 (ln/i) 5 / 3 . (i 

(102) a, ~ , 77* ~ ^— ^ — , A* - 
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At the transition point A t (/i), the equality of the entropy of the two phases leads to the 
relation 

3 3 
(103) - - In rj cond - 3a In // - (1 - a) ln(l - a) = --hi7? ffas . 

In the limit /t — > +00, the fraction of mass contained in the condensate behaves like 

1 



(104) a ~ 
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The energy of transition and the temperature of each phase are given by 

/a (ln/i) 5 / 3 
(105) A t — %as ~ r/ cond . 

The jump of temperature at the transition point is 

1 1 n 



(106) 



Vcond Vgas (ln/i) 5 / 3 ' 



For a given energy, the thermodynamical parameters of the condensate behave, in the limit 
\i — > +00, like 

fin* 1 (ln/i) 5/3 q . 

m/i // 

For the germ, we have 

(108) CXgerm ~ ( ~ J > Vgerrn * Vgas: ^ germ * $g as . 

In the canonical ensemble, the maximum temperature of the condensed phase, corresponding 
to 77* (//) , is determined by the relation 

(109) \ma- l l* - — ?— - ^ = 0. 
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In the limit fx — > +00, the thermodyanical parameters behave at the point of maximum tem- 
perature like 

1 In /i /i 
(110) l-ct*~- , 77* ~ , A* - A maa .(/x) 



m /x /i m fx 

At the transition point i] t (fx), the equality of the free energy of the two phases yields 

3 3 3 3 

(111) -n = -3a In w - (1 - a) ln(l - a) + -fir]a 5/3 + -nail - a) + -7/(1 - a) 2 . 
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In the limit fx — ► +00, the temperature of transition behaves like 

In fx 
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and the fraction of mass contained in the nucleus like 

(113) i_ a „J_ 

fx 2 

The energy of the two phases at the transition point is given by 

(114) A gas ~ — -^— , A cond — A max (ji) 



In /i in /i 

and the jump of energy is 

(115) A cond - A gas ~ /i. 

For a given temperature, the thermodynamical parameters characterizing the condensate 
behave, in the limit fx — > +00, like 

(116) 1 - a cond ~ e~ m , A ~ A maa .(/i) ~ fx, J ~ i]A max (fx) ~ //. 
For the germ, we have 

1 \ 3 / 2 
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(117) ftgerm ~ ^"/^J ' ^ ~^ ^ SaS ' ^ ~^ ^ 9aS ' 

6.4 Comparision with previous works 

Phase transitions in self-gravitating systems were first investigated with the aid of toy models 
which could be solved exactly without recourse to a meanfield approximation. For example, 
Lynden-Bell & Lynden-Bell |13[ considered a system of N particles confined to the surface of 



a sphere of variable radius. They calculated exactly the density of states in the microcanonical 
ensemble and showed the existence of a region with negative specific heats. Then, they evaluated 
the partition function in the canonical ensemble and demonstrated that the region of negative 
specific heats is replaced by a remarkable giant phase transition connecting a "gaseous" phase 
(at high energies) to a "condensed" phase (at low energies). Padmanabhan |I[ obtained similar 
results with a simpler model consisting on only two particles in gravitational interaction confined 
within a spherical box. The phase diagram determined by these authors is similar to the one 



reported in Fig. 10. These models exhibit a phase transition in the canonical ensemble but not 
in the microcanonical ensemble (unlike in Fig. |7|). This is because, for these simple models, 
the density of states remains finite when the small-scale cut-off a is set equal to zero whereas 
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for more realistic self-gravitating systems, it diverges. By contrast, the partition function 
is divergent for a = and this leads to the occurence of a phase transition in the canonical 
ensemble when a is sufficiently small. Padmanabhan investigated the dependance of the critical 
parameters with the small scale cut-off a. In particular, he found that the temperature of 
transition is given by 

^ Cm 2 
118 T t 



3a \n(R/a) 



This expression qualitatively agrees with our result ( |1 1 2j ) , which becomes in dimensional vari- 
ables 
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Recall that N = 2 in Padmanabhan's model. He also computed the change of energy at the 
transition point and found that 

Gm 2 ( 1 
(120) E gas — E cond = 1 1 — 



3ki(R/a) 

This expression also qualitatively agrees with our result ( |115| ), leading to 

GN 5/3 m 2 

(121) E gas — E cond ~ . 

a 

The logarithmic correction in Eq. ( |1 20|) is a particularity of Padmanabhan's model arising 
from the low value of N . Finally, Padmanabhan investigated the dependance of the T(E) 
curve with the small-scale cut-off a (see his Fig. 3.2). For large a, his diagram is similar 
to that of Fig. [TT]. In particular, there exists a critical short distance cut-off above which 
the phase transition disappears (tricritical point). It should be stressed that the statistical 
approach based on the evaluation of g(E) or Z{(3) does not determine the metastable states 
unlike the thermodynamical approach based on the maximization of S[p] or J[p\. Only the true 
equilibrium states, which correspond to global maxima of S or J, appear in the T(E) diagram 
(in other words, we directly obtain the "plateaux" without beeing required to make a Maxwell 
construction). These equilibrium states are expected to be reached for t — ► +oo but they are 
not necessarily the most relevant for astrophysical applications: as discussed previously, the 
metastable equilibrium states may be long-lived and may correspond to the structures that 
are actually observed in the universe. Indeed, the statistical mechanical approach tells nothing 
about the timescales involved in the establishement of the equilibrium. 

Phase transitions in self-gravitating systems have also been investigated in the meanfield 
approximation for less idealized models. Aronson & Hansen M have considered the case of a 
classical hard spheres gas modeled by a van der Waals equation of state. They considered a 
relatively large cut-off and evidenced only the phase transition in the canonical ensemble (they 
obtained a diagram similar to that of Fig. |10|) . They also proposed a simple analytical model 
to describe this phase transition. In their model, the "gaseous" phase consists of iV particles 
spread with uniform density throughout the whole container and the "condensed" phase has 
all N particles collapsed into a central core of uniform density. The model that we proposed in 
Sec. |6] is more general because we allow the condensate to contain an arbitrary fraction of the 
total mass. Then, the fraction of mass that is actually achieved at equilibrium is determined by 
maximizing the free energy versus a. In the canonical ensemble (the only situation discussed 
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by Aronson & Hansen), the fraction of mass contained in the condensate is close to one, so 
that our precedure provides additional support to their Ansatz. However, our model allows 
us to describe also the unstable solution with the "germ" (a < 1) and to obtain a better 
representation of the whole bifurcation diagram (see Fig. |2"3] ). In addition, our model can also 
describe the phase transitions in the microcanonical ensemble (occuring for sufficiently small 
values of the cut-off a) which was not considered by Aronson and Hansen. 
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Figure 23: Equilibrium phase diagram for a hard spheres gas obtained in the framework of our 
analytical model. The parameters have been chosen so as to correspond to Fig. 2 of Aronson & 
Hansen [H]: R = 60 km, R = 30 km, R = 15 km; N = 10 57 ; m = neutron mass = 1.67 x 10~ 24 g; 
a = 0.4 x 10~ 13 cm. These parameters correspond to a neutron-star-like structure. In our 
analytical model, we have adopted a value of the geometrical factor g = 2 The agreement 
with the numerical study of Aronson & Hansen is fairly good. For small radius R = 15 km, 
there is no phase transition. For larger radii (R = 30 km, R = 60 km), a phase transition 
connects the "gaseous" phase (upper branch) to the "condensed" phase (lower branch). This 
phase transition forms a Maxwell plateau (dashed-line) at a temperature depending on the size 
R of the system (according to Eqs. ( |111| ) and (^)). In dot-dashed line, we have ploted the 
free energy F = E — TS as a function of the inverse temperature for the case R = 60 km. The 
crossing point determines the temperature of transition. 

The study of Aronson & Hansen was reconsidered by Stahl et al. || who pointed out that 
the van der Waals equation of state does not adequately describe hard spheres systems at high 
densities. They considered a more general equation of state and studied in detail the case of 
small filling factors for which a phase transition occurs in the microcanonical ensemble. They 
also determined numerically the dependance of the transition temperature (in the canonical 
ensemble) as a function of the filling factor and considered applications of their results in the 
context of planet formation. 

Following the works of Aronson and Hansen || and Stahl et al. ||, different authors have 
attempted to describe phase transitions in self-gravitating systems by introducing a small scale 
regularization of the gravitational potential. For example, in the study of Follana & Laliena 
1T0| , the softening is achieved by truncating to N terms an expansion of the Newtonian potential 



in spherical Bessel functions. These authors obtained an equilibrium phase diagram similar to 
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the one in Fig. [K]. However, they could not achieve very large values of N in their study so that 
they were unable to see the development of the spiral (and the corresponding phase transition 
in the microcanonical ensemble) as iV — > +00. Another analysis of phase transitions in self- 
gravitating systems was provided by Youngkins & Miller [O] with a one- dimensional model of 
concentric spherical mass shells. They studied this system in the microcanonical, canonical and 
grand canonical ensembles both numerically and analytically in the meanfield approximation. 
They found an overall good agreement between their numerical simulations and the meanfield 
predictions. They observed phase transitions in the microcanonical and canonical ensembles 
but not in the grand canonical ensemble in which the system remains homogeneous. This last 
result may be, however, an artefact of their one-dimensional model. It is plausible that, in the 
grand canonical ensemble, the self-gravitating gas fragments in a series of clumps (at different 
scales) as observed in cosmology and for the interstellar medium. Some theoretical arguments 



in favour of this scenario have been given by Semelin et al. |30j and Chavanis |25| . Of course, to 
study the development of these clumps, it is necessary to extend the thermodynamical analysis 
to the full three dimensional problem and relax the assumption of spherical symmetry. 

The thermodynamics of self-gravitating fermions was investigated by Lynden-Bell & Wood 
and Hertel & Thirring in the early seventies (see the discussion of Aronson & Hansen ||), 
but these papers were apparently not published. In Refs. ||, |7|], it is proved that a rigorous 
thermodynamic limit exists for self-gravitating fermions but the corresponding phase diagram 
is not explicitly given. This equilibrium phase diagram was calculated by Bilic & Viollier |J for 
a particular value of the degeneracy parameter adapted to a cosmological settling. It was also 
calculated independantly by Chavanis & Sommeria || for an arbitrary degree of degeneracy in 
the context of the theory of "violent relaxation" . In Ref . , the development of the spiral for 
high values of the degeneracy parameter and the associated phase transition that occurs in the 
microcanonical ensemble (in addition to the more well-known phase transition in the canonical 
ensemble) was clearly shown. However, the point of transition was not explicitly determined 
and this has been done in the present study. 



6.5 Astrophysical applications 



The application of the hard spheres model in astrophysics could concern the fragmentation 
of the interstellar medium and the formation of stars or even smaller objects such as planets 
||. These objects would correspond to the "condensate" that results from the collapse of a 
cloud of gas or dust. On the other hand, the model of self-gravitating fermions could have 
applications for massive neutrinos in dark matter models [|16|, |5| , white dwarfs |14| and neutron 
stars |[5L M. By cooling below a critical temperature, a condensed phase emerges consisting 



of a completely degenerate nucleus surrounded by a dilute enveloppe, as extensively studied 
in early models of stellar structure ||14| . This model could also be relevant for the "violent 
relaxation" of collisionless stellar systems [] such as elliptical galaxies [17, 20]. In that 
case, the exclusion principle is a consequence of the Liouville theorem. Since degeneracy can 
stabilize the system without changing its overall structure at large distances, we have suggested 
in Ref. || that degeneracy could play a role in galactic nuclei. The recent simulations of 
Leeuwin & Athanassoula 
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and the theoretical model of Stiavelli |33| are consistent with 
this idea especially if the nucleus of elliptical galaxies contains a primordial massive black 
hole. Indeed, the effect of degeneracy (in the sense of Lynden-Bell) on the distribution of stars 
surrounding the black hole can explain the cusps observed at the center of galaxies. Whether 



Hn fact, the problem is complicated because violent relaxation eventually fades before the maximum entropy 
state is attained. Thus, equation (|J) [or Eq. (0)] is unlikely to be reached throughout the whole cluster. However, 
it is reasonable to hold in the central region in which violent relaxation occurs most violently jl?], . 
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or not elliptical galaxies are degenerate remains however a matter of debate because when the 
core becomes dense, two-body encounters will come into play and break the Liouville theorem 
(Shu's criticism ||34j| ). This form of degeneracy may, however, be relevant for massive neutrinos 
in Dark Matter models where it competes with quantum degeneracy [|TjJ. In fact, the thermal 
equilibrium distribution of massive neutrinos in Dark Matter models might be justified more 
by the process of "violent relaxation" than by a collisional relaxation. Indeed, the time scale of 
gravitational two-body encounters for neutrinos is extremely long so that the criticism raised 
by Shu does not apply |35, |I9 |. Therefore, the commonly adopted Fermi-Dirac distribution 
of self-gravitating neutrinos might be due to Lynden-Bell's type of degeneracy rather than to 
quantum mechanics. Anyway, whatever the source of the exclusion principle (Lynden-Bell or 
Pauli), the self-gravitating Fermi-Dirac model predicts the formation of a dense degenerate 
nucleus ("fermion star") with a small radius and a large mass ||, ||. As suggested in Refs. 
|| H, |36||, this dense degenerate nucleus could be an alternative to black holes at the center of 
galaxies. On the other hand, at large distances, the density of the (isothermal) self-gravitating 
Fermi gas decays like r~ 2 which is a condition that dark galactic halos must fulfill in order to 



reproduce the flat rotation curves of spiral galaxies |[22|| . Therefore, this model of self-gravitating 



fermions has a chance to account for the structure of dark matter in galactic halos. 



7 Conclusion 

In this paper, we have described the inequivalence of statistical ensembles and the nature of 
phase transitions in self-gravitating systems by considering the case of self-gravitating fermions 
or the case of a classical hard spheres gas. The introduction of an effective repulsion at short 
distances avoids the singularity of the "naked" gravitational potential. It is likely that similar 
results will be obtained with different forms of the regularization. For large values of the cut-off 
a, there are no phase transition. For intermediate values of a, phase transtions occur in the 
canonical ensemble but not in the microcanonical ensemble. The corresponding phase diagram 
is of the type of Fig. [LT] and has been found by various authors || |5|, || [LTJ . For smaller values 
of a, phase transitions occur both in the canonical and in the microcanonical ensemble. The 
corresponding phase diagram is of the type of Fig. [7] and was first obtained in Ref. [0. As the 
small-scale cut-off is decreased, the T(E) curve winds up and tends to the classical spiral for 
a -> (see Fig. 0). 

Depending on the value of the cut-off a and of the ensemble considered (microcanonical or 
canonical), three kinds of phase transitions can be evidenced which separate a "gaseous" phase 
from a "condensed" phase. In the microcanonical ensemble and for sufficiently small a (Fig. 
0), a gravitational first order phase transition occurs at an energy E t (a) at which the "gaseous" 
states pass from global to local entropy maxima and the "condensed" states from local to global 
entropy maxima. This transition is marked by a discontinuity of the temperature and of the 
specific heats. In the limit a — > 0, the transition energy is rejected to +oo, so that the stable 
(gaseous) solutions of the classical spiral are only metastable. These states may be physical 
depending whether the initial condition lies in their "basin of attraction" or not || [0], ^6|| . 
However, the metastable branch disappears at a critical energy E c , discovered by Antonov |2j, 
at which the "gravothermal catastrophe" || occurs. This collapse can be considered as a zeroth 
order phase transition since it is associated with a discontinuity of entropy. For larger values 



of a, there is only one entropy maximum for each value of energy (Fig 10) and the previously 
described phase transitions are suppressed. However, as we progressively decrease energy, the 
self-gravitating gas acheives higher and higher density contrasts and builds up a compact core 
containing more and more mass. This gravitational "clustering" can be called a second order 
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phase transition since the specific heats diverges at the critical point and the order parameter 
experiences a rapid variation in the region of negative specific heats (but remains continuous). 
Similar phase transitions occur in the canonical ensemble. It is interesting that these phase 
transitions can be understood with the aid of a simple analytical model which allows one to 
determine the dependance of the thermodynamical parameters with the cut-off value. The 
present study can be extended to include rotation (in preparation). It is also important to 
develop non equilibrium models to determine the structure of the "basin of attraction" of 
the equilibrium states when several solutions exist. A first step in that direction was made by 
Youngkins & Miller [TT] and by Chavanis et al. |^B1 with the aid of simplified dynamical models. 
These models could be used in particular to investigate numerically the occurence of the phase 
transitions and the robustness of the metastable states. If the system is placed in a metastable 
state (local entropy maximum), it will eventually jump to the global entropy maximum, but 
this can take an infinite (physically irrelevant) time. Indeed, the probability that a fluctuation 
will allow the phase transition to develop is expected to be extremely low (except maybe near 
the critical point) |37| . Therefore, if the system is trapped in a metastable state (but still slowly 
evolving along the series of equilibrium by losing mass or energy like for globular clusters), the 
phase transition will occur at the critical point A c rather than at A t . However, if the system 
is initially far from equilibrium, there is no simple criterion to decide a priori whether it will 
converge towards the local or the global entropy maximum. Only direct numerical simulations 
can answer this question and sketch the structure of the basin of attraction for self-gravitating 
systems. 
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A Entropy of the self-gravitating Fermi gas 



In this Appendix, we give the main steps for deriving the expression ( |L7D for the entropy. 
Substituting the Fermi-Dirac distribution 

(122) / ,/n 



v2_ ' 



1 + ke^e^- 

in the entropy (|5|), we obtain after some rearrangements 

(123) r] S = M\nk + J p^dh + Kfi + rjo J\n(l + ^e^e- p! £)d 3 rd s v. 
The last integral can be integrating by parts yielding the value 20K/3. Therefore, 

(124) r] S = M In k + J pipdh + 
Using the definition of ip, we get 

(125) r] S = M\nk + 2(3W - M(3% + -K(3. 

3 

Now, the central density is determined by the relation ip{a) = /3($(R) — $o) with $(i?) = 
-GM/R. Hence 

(126) i* =]uk + ?l W+ WL + H a)+ *JX. 
K 1 M M R ^ K ' 3M 

Using Eq. (3.12) of Chavanis & Sommeria to express the potential energy W and the kinetic 
energy K = E — W in terms of E, we find 
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(127) ^ = \nk + V + tlj(a) - -At? - -r]^-I 3/2 (ke^ a) ). 



M y ' 3 ' 9 ' fi 

Finally, using r] = fi 2 /a 4 , resulting from Eqs. (|i~3"D and (|l^), we obtain Eq. (p!7|). 
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